##### This script recreates the article figures and tables reported in the online appendix.


#####################################################################################################
## Figure 1 (Vývoj průměrné podpory přerozdělování příjmů a míry ekonomické aktivity v 27 zemích.) ##
#####################################################################################################

### Aggregating the data for the country-round time series figure
ess_cntry_round_figure <- ess_4_ld_models[, .(wght_mean_support = weighted.mean(x = red_sup, w = pspwght),
                                              LFPR = mean(LFPR_y0), dev_mean_LFPR_y0 = mean(dev_mean_LFPR_y0)), 
                                          by = c("cntry", "essround")]

countries <- unique(ess_cntry_round_figure$cntry)

### a. COUNTRY-ROUND SUPPORT FOR REDISTRIBUTION
# creating an empty matrix for the purposes of creating the time series plots
full_matrix_red_sup <- matrix(data = rep(NA, 27*9), nrow = 27, ncol = 9)
rownames(full_matrix_red_sup) <- countries
colnames(full_matrix_red_sup) <- c("Round 1", "Round 2", "Round 3", "Round 4", "Round 5", 
                                   "Round 6", "Round 7", "Round 8", "Round 9")

# assignment of available aggregate country-year support for redistribution
for(i in 1:nrow(ess_cntry_round_figure)){
  country_i <- ess_cntry_round_figure$cntry[i]
  round_i <- as.numeric(ess_cntry_round_figure$essround[i])
  
  row_full <- as.numeric(which(rownames(full_matrix_red_sup) == country_i))
  column_full <- round_i
  
  full_matrix_red_sup[row_full, column_full] <- ess_cntry_round_figure$wght_mean_support[i]
}

### b. COUNTRY-ROUND LABOUR FORCE PARTICIPATION RATE
# creating an empty matrix for the purposes of creating the time series plots
full_matrix_LFPR <- matrix(data = rep(NA, 27*9), nrow = 27, ncol = 9)
rownames(full_matrix_LFPR) <- countries
colnames(full_matrix_LFPR) <- c("Round 1", "Round 2", "Round 3", "Round 4", "Round 5", 
                                "Round 6", "Round 7", "Round 8", "Round 9")

# assignment of available aggregate country-year support for redistribution
for(i in 1:nrow(ess_cntry_round_figure)){
  country_i <- ess_cntry_round_figure$cntry[i]
  round_i <- as.numeric(ess_cntry_round_figure$essround[i])
  
  row_full <- as.numeric(which(rownames(full_matrix_LFPR) == country_i))
  column_full <- round_i
  
  full_matrix_LFPR[row_full, column_full] <- ess_cntry_round_figure$LFPR[i]
}

# range of average support for income redistribution across 202 country rounds
range(ess_cntry_round_figure$wght_mean_support)

# range of LFPR across 202 country rounds
range(ess_cntry_round_figure$LFPR)

full_names <- c("Rakousko (AT)", "Belgie (BE)", "Bulharsko (BG)", "Švýcarsko (CH)", "Kypr (CY)", 
                "Česká republika (CZ)", "Německo (DE)", "Dánsko (DK)", "Estonsko (EE)", "Španělsko (ES)", 
                "Finsko (FI)", "Francie (FR)", "Spojené království (GB)", "Řecko (GR)", "Maďarsko (HU)", 
                "Irsko (IE)", "Izrael (IL)", "Island (IS)", "Itálie (IT)", "Litva (LT)", 
                "Nizozemsko (NL)", "Norsko (NO)", "Polsko (PL)", "Portugalsko (PT)", "Švédsko (SE)", 
                "Slovinsko (SI)", "Slovensko (SK)")

##### A Creating the time series plot: countries 1:15 ranked alphabetically based on country code
### Exporting the plot: (PDF Size = Device Size: 12 x 8)
par(mar = c(2.5, 2.5, 2.5, 2.5))
mat <- matrix(c(1:15),
              nrow = 3, ncol = 5,
              byrow = TRUE)
layout(mat = mat)

# Both country names and codes
for(i in 1:15){
  plot(seq(from = 1, to = 9, by = 1), full_matrix_red_sup[i, ], ylim = c(1.75, 3.9), main = full_names[i], 
       xaxt = "n", pch = 19, cex = 1.3, type = "o", xlab = "", ylab = "", lty = 3, lwd = 1.5, las = 1)
  axis(1, at = 1:9, labels = c(2002,2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018))
  correlation <- round((cor(full_matrix_red_sup[i, ], full_matrix_LFPR[i, ], use = "complete.obs")), digits = 2)
  text(x = 2, y = 3.775, labels = bquote(r == .(correlation)), cex = 1.1)
  par(new = T)
  plot(seq(from = 1, to = 9, by = 1), full_matrix_LFPR[i, ], ylim = c(55, 95), xaxt = "n", yaxt="n", pch = 15, 
       cex = 1.3, type = "o", xlab = "", ylab = "", lty = 1, col = "gray")
  axis(side = 4, las = 1, col.axis = "gray", col.ticks = "gray")
}

##### B Creating the time series plot: 16:27 countries ranked alphabetically based on country code
### Exporting the plot: (PDF Size = Device Size: 12 x 8)
mat <- matrix(c(1:12, 13, 13, 13),
              nrow = 3, ncol = 5,
              byrow = TRUE)
layout(mat = mat)

# Both country names and codes
for(i in 16:27){
  plot(seq(from = 1, to = 9, by = 1), full_matrix_red_sup[i, ], ylim = c(1.75, 3.9), main = full_names[i], 
       xaxt = "n", pch = 19, cex = 1.3, type = "o", xlab = "", ylab = "", lty = 3, lwd = 1.5, las = 1)
  axis(1, at = 1:9, labels = c(2002,2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018))
  correlation <- round((cor(full_matrix_red_sup[i, ], full_matrix_LFPR[i, ], use = "complete.obs")), digits = 2)
  text(x = 2, y = 3.775, labels = bquote(r == .(correlation)), cex = 1.1)
  par(new = T)
  plot(seq(from = 1, to = 9, by = 1), full_matrix_LFPR[i, ], ylim = c(55, 95), xaxt = "n", yaxt="n", pch = 15, 
       cex = 1.3, type = "o", xlab = "", ylab = "", lty = 1, col = "gray")
  axis(side = 4, las = 1, col.axis = "gray", col.ticks = "gray")
}

plot.new()
legend("top", pch = c(19, 15), col = c("black", "gray"), lty = c(3, 1), 
       legend = c("Podpora přerozdělování příjmů (vertikální osa vlevo)", "Míra ekonomické aktivity (v %, vertikální osa vpravo)"), 
       bty = "n", cex = 1.75, pt.cex = c(1.75, 1.75), text.col = c("black", "gray"), y.intersp = 1.75)



###########################################################################################################
## Figure 2 (Longitudinální a průřezový vztah mezi podporou přerozdělování a mírou ekonomické aktivity.) ##
###########################################################################################################

# Aggregated data for the COUNTRY-ROUND effect figure (WE)
par(mfrow = c(2, 1), mar = c(4.25, 4.5, 2.25, 1))
# Exporting the plot: (PDF Size = Device Size: 9.00 x 7.95)

# Plot to the article
plot(ess_cntry_round_figure$dev_mean_LFPR_y0, ess_cntry_round_figure$wght_mean_support, pch = 19, ylim = c(1.75, 4), xlim = c(-10, 10), 
     xlab = expression(paste("Odchylka od průměrné míry ekonomické aktivity v zemi", ~ (LFPR[zrCWC]))), 
     ylab = expression(paste("Podpora přerozdělování (země-rok)", ~~~ (bar(y)[zr]))), cex.lab = 0.85,
     las = 1, main = "Longitudinální vztah")
abline(lm(ess_cntry_round_figure$wght_mean_support ~ ess_cntry_round_figure$dev_mean_LFPR_y0), lwd = 2.5, lty = 2, col = "grey")
cor(ess_cntry_round_figure$dev_mean_LFPR_y0, ess_cntry_round_figure$wght_mean_support, use = "complete.obs")
text(x = 8, y = 1.9, labels = expression(paste("r = 0.01, ", ~ n[zr] == 202)), cex = 1.25)
mtext(expression(bar(y)[zr] == 2.857 + 0.002*LFPR[zrCWC]), side = 3, line = -2,  col = "gray", cex = 1.1)
lm(ess_cntry_round_figure$wght_mean_support ~ ess_cntry_round_figure$dev_mean_LFPR_y0)

# Aggregated data for the COUNTRY effects figures (BE)
# computing the weighted mean-country support for redistribution ("mean_support_w") - each weighted mean is based on 4-9 country-round means (i.e. this procedure aggregates the country-round means)
ess_cntry_figure <- ess_cntry_round_figure[, .(mean_support_w = mean(wght_mean_support), mean_LFPR = mean(LFPR)), by = .(cntry)]

# Plot to the article
plot(ess_cntry_figure$mean_LFPR, ess_cntry_figure$mean_support_w, pch = 19, ylim = c(1.8, 4), xlim = c(60, 90), 
     xlab = expression(paste("Průměrná míra ekonomické aktivity v zemi", ~~ (bar(LFPR)[z]))), 
     ylab = expression(paste("Průměrná podpora přerozdělování v zemi ", ~~ (bar(y)[z]))), las = 1, cex.lab = 0.85,
     main = "Průřezový vztah", cex = 1.2)
abline(lm(ess_cntry_figure$mean_support_w ~ ess_cntry_figure$mean_LFPR), lwd = 2, lty = 2, col = "gray")
text(ess_cntry_figure$mean_LFPR, ess_cntry_figure$mean_support_w, labels = ess_cntry_figure$cntry, pos = 3, cex = 0.7)
cor(ess_cntry_figure$mean_LFPR, ess_cntry_figure$mean_support_w)
text(87, 2, labels = expression(paste("r = -0.53, ", ~ n[z] == 27)), cex = 1.25)
mtext(expression(bar(y)[z] == 5.165 -0.031*bar(LFPR)[z]), side = 3, line = -2.25,  col = "gray", cex = 1.1)
lm(ess_cntry_figure$mean_support_w ~ ess_cntry_figure$mean_LFPR)



#######################################################################################################
## Figure 3 (Graf 3. Bodové odhady a 95% intervaly spolehlivosti pro efekty míry ekonomické aktivity ##
## při různých konfiguracích víceúrovňových modelů.) ##
#######################################################################################################

par(mar = c(4.2, 0.5, 4.2, 0.5), mfcol = c(3, 2))
# PDF Size = Device size: 12,5 x 8 (Portrait)

# Figure i: models controlling for compositional effects
# All explanatory variables entered in their original form
plot(x = c(-0.06, 0.0425), y = c(7.2, 0.8), type = "n", ylab = "", xlab = "",  yaxt = "n", xaxt = "n")
title(main = "i: Kompoziční efekty (Tabulka OA1)", line = 2.75, cex.main = 1.25)
title(xlab = "Koeficient: míra ekonomické aktivity", line = 2.25, cex.lab = 1.1)

points(x = LFPR_initial$Estimate, y = 7:1, pch = 19, cex = 1)

for(i in 1:nrow(LFPR_initial)){
  lines(x = c(LFPR_initial$lower_bound[i], LFPR_initial$upper_bound[i]), y = c(8-i, 8-i), lwd = 1.5)
  text(x = -0.05, y = 8-i, labels = paste(c(LFPR_initial[i, 1], " (", round(LFPR_initial[i, 4], digits = 1), ")"), collapse = ""))
}
abline(v = 0, lty = 2)
axis(side = 3, at = seq(from = -0.04, to = 0.02, by = 0.01))
axis(side = 1, at = seq(from = -0.04, to = 0.02, by = 0.01))


# Figure j: models NOT controlling for compositional effects
# Individual-level explanatory variables centered within country rounds (CWC)
plot(x = c(-0.06, 0.0425), y = c(7.2, 0.8), type = "n", ylab = "", xlab = "",  yaxt = "n", xaxt = "n")
title(main = "j: Bez kompozičních efektů  (Tabulka OA2)", line = 2.75, cex.main = 1.25)
title(xlab = "Koeficient: míra ekonomické aktivity", line = 2.25, cex.lab = 1.1)

points(x = LFPR_CWC_individual$Estimate, y = 7:1, pch = 19, cex = 1)

for(i in 1:nrow(LFPR_CWC_individual)){
  lines(x = c(LFPR_CWC_individual$lower_bound[i], LFPR_CWC_individual$upper_bound[i]), y = c(8-i, 8-i), lwd = 1.5)
  text(x = -0.05, y = 8-i, labels = paste(c(LFPR_CWC_individual[i, 1], " (", round(LFPR_CWC_individual[i, 4], digits = 1), ")"), collapse = ""))
}
abline(v = 0, lty = 2)
axis(side = 3, at = seq(from = -0.04, to = 0.02, by = 0.01))
axis(side = 1, at = seq(from = -0.04, to = 0.02, by = 0.01))

# Empty plot
plot.new()


# Figure k: Models controlling for compositional effects, but NOT controlling for other contextual variables
plot(x = c(-0.0725, 0.055), y = c(7.2, 0.6), type = "n", ylab = "", xlab = "",  yaxt = "n", xaxt = "n")
title(main = "k: Kompoziční efekty: efekt LFPR mezi zeměmi (BCE) a uvnitř zemí (WCE) (Tabulka OA3)", line = 2.75, cex.main = 1.25)
title(xlab = "Koeficient: míra ekonomické aktivity", line = 2.25, cex.lab = 1.1)

points(x = only_LFPR_BE$Estimate, y = (7:1)+0.1, pch = 19, cex = 1)

for(i in 1:nrow(only_LFPR_BE)){
  lines(x = c(only_LFPR_BE$lower_bound[i], only_LFPR_BE$upper_bound[i]), y = c(8-i, 8-i)+0.1, lwd = 1.5)
  text(x = -0.054, y = 8-i, labels = paste(c("BCE ", only_LFPR_BE[i, 1], " (", round(only_LFPR_BE[i, 4], digits = 1), ")"), collapse = ""))
}
abline(v = 0, lty = 2)
axis(side = 3, at = seq(from = -0.04, to = 0.02, by = 0.01))
axis(side = 1, at = seq(from = -0.04, to = 0.02, by = 0.01))

points(x = only_LFPR_WE$Estimate, y = (7:1)-0.25, pch = 19, cex = 1, col = "gray60")

for(i in 1:nrow(only_LFPR_WE)){
  lines(x = c(only_LFPR_WE$lower_bound[i], only_LFPR_WE$upper_bound[i]), y = c(8-i, 8-i)-0.25, lwd = 1.5, col = "gray60")
  text(x = 0.045, y = (8-i)-0.225, labels = paste(c("WCE ", only_LFPR_WE[i, 1], " (", round(only_LFPR_WE[i, 4], digits = 1), ")"), collapse = ""), col = "gray60")
}


# Figure l: models NOT controlling for compositional effects, but splitting the contextual effect into BE and WE
# Individual-level explanatory variables centered within country rounds (CWC)
plot(x = c(-0.0725, 0.055), y = c(7.2, 0.6), type = "n", ylab = "", xlab = "",  yaxt = "n", xaxt = "n")
title(main = "l: Bez kompozičních efektů: efekty mezi zeměmi (BCE) a uvnitř zemí (WCE) (Tabulka OA4)", line = 2.75, cex.main = 1.25)
title(xlab = "Koeficient: míra ekonomické aktivity", line = 2.25, cex.lab = 1.1)

points(x = LFPR_CWC_ind_BE$Estimate, y = (7:1)+0.1, pch = 19, cex = 1)

for(i in 1:nrow(LFPR_CWC_ind_BE)){
  lines(x = c(LFPR_CWC_ind_BE$lower_bound[i], LFPR_CWC_ind_BE$upper_bound[i]), y = c(8-i, 8-i)+0.1, lwd = 1.5)
  text(x = -0.059, y = 8-i, labels = paste(c("BCE ", LFPR_CWC_ind_BE[i, 1], " (", round(LFPR_CWC_ind_BE[i, 4], digits = 1), ")"), collapse = ""))
}
abline(v = 0, lty = 2)
axis(side = 3, at = seq(from = -0.04, to = 0.02, by = 0.01))
axis(side = 1, at = seq(from = -0.04, to = 0.02, by = 0.01))

points(x = LFPR_CWC_ind_WE$Estimate, y = (7:1)-0.25, pch = 19, cex = 1, col = "gray60")

for(i in 1:nrow(LFPR_CWC_ind_WE)){
  lines(x = c(LFPR_CWC_ind_WE$lower_bound[i], LFPR_CWC_ind_WE$upper_bound[i]), y = c(8-i, 8-i)-0.25, lwd = 1.5, col = "gray60")
  text(x = 0.036, y = (8-i)-0.225, labels = paste(c("WCE ", LFPR_CWC_ind_WE[i, 1], " (", round(LFPR_CWC_ind_WE[i, 4], digits = 1), ")"), collapse = ""), col = "gray60")
}

# Figure m: models controlling for compositional effects and splitting the contextual effect into BE and WE
plot(x = c(-0.0725, 0.055), y = c(7.2, 0.6), type = "n", ylab = "", xlab = "",  yaxt = "n", xaxt = "n")
title(main = "m: Kompoziční efekty: efekty mezi zeměmi (BCE) a uvnitř zemí (WCE) (Tabulka OA5)", line = 2.75, cex.main = 1.25)
title(xlab = "Koeficient: míra ekonomické aktivity", line = 2.25, cex.lab = 1.1)

points(x = LFPR_BE$Estimate, y = (7:1)+0.1, pch = 19, cex = 1)

for(i in 1:nrow(LFPR_BE)){
  lines(x = c(LFPR_BE$lower_bound[i], LFPR_BE$upper_bound[i]), y = c(8-i, 8-i)+0.1, lwd = 1.5)
  text(x = -0.054, y = 8-i, labels = paste(c("BCE ", LFPR_BE[i, 1], " (", round(LFPR_BE[i, 4], digits = 1), ")"), collapse = ""))
}

abline(v = 0, lty = 2)
axis(side = 3, at = seq(from = -0.04, to = 0.02, by = 0.01))
axis(side = 1, at = seq(from = -0.04, to = 0.02, by = 0.01))

points(x = LFPR_WE$Estimate, y = (7:1)-0.25, pch = 19, cex = 1, col = "gray60")

for(i in 1:nrow(LFPR_WE)){
  lines(x = c(LFPR_WE$lower_bound[i], LFPR_WE$upper_bound[i]), y = c(8-i, 8-i)-0.25, lwd = 1.5, col = "gray60")
  text(x = 0.04, y = (8-i)-0.225, labels = paste(c("WCE ", LFPR_WE[i, 1], " (", round(LFPR_WE[i, 4], digits = 1), ")"), collapse = ""), col = "gray60")
}

# Gray rectangle: Model D
rect(xleft = -0.0775, ybottom = 3.355, xright = 0.061, ytop = 4.355, 
     col = rgb(0.3, 0.3, 0.3, alpha = 0.125), border = "black")

# Gray rectangle: Model F
rect(xleft = -0.0775, ybottom = 1.355, xright = 0.061, ytop = 2.355, 
     col = rgb(0.3, 0.3, 0.3, alpha = 0.125), border = "black")



#############################
#####  ONLINE APPENDIX  #####
#############################

library(data.table)

##### Table A1: Actual sample sizes
table_A1 <- table(ess_4_ld_models$cntry, ess_4_ld_models$essround)
dim(ess_4_ld_models)
write.csv(x = table_A1, file = "Table_A1_appendix.csv", dec = ".")

# Corresponding table when ignoring missing values
table_A1_full <- table(ess_4$cntry, ess_4$essround)
table_A1_full["FR", "1"] <- 0
table_A1_full["FR", "2"] <- 0
nrow(ess_4)
sum(table_A1_full)

# Proportion of excluded respondens due to missing values
sum(table_A1)/sum(table_A1_full)


##### Table A2: Country means of individual-level explanatory variables
describe_individual <- ess_4_ld_models[, .(employed = mean(employed), unemployed = mean(unemployed), hinc_2_cop = mean(hinc_2_cop), hinc_3_dif = mean(hinc_3_dif), hinc_4_vdif = mean(hinc_4_vdif),
                                           female = mean(female), age = mean(age), isced2 = mean(isced2), isced3 = mean(isced3), isced4 = mean(isced4), isced56 = mean(isced56), stfeco = mean(stfeco), lrscale = mean(lrscale)), 
                                       by = "cntry"]

write.csv(x = describe_individual, file = "Table_A2_appendix.csv", dec = ".")

# Country-round means of individual-level explanatory variables (not reported in the online appendix of the article)
ess_4_ld_models[, .(employed = mean(employed), unemployed = mean(unemployed), hinc_2_cop = mean(hinc_2_cop), hinc_3_dif = mean(hinc_3_dif), hinc_4_vdif = mean(hinc_4_vdif),
                    female = mean(female), age = mean(age), isced2 = mean(isced2), isced3 = mean(isced3), isced4 = mean(isced4), isced56 = mean(isced56), stfeco = mean(stfeco), lrscale = mean(lrscale)), 
                by = "cntry_round"]


##### TABLE A3: Values of contextual variables at country-round level
table_A3 <- ess_4_ld_models[, .(quarter = key_quartal[1], LFPR = mean(LFPR_y0), SOC_EXP = mean(SOCEXP_y0), GDP_capita = mean(GDP_y0), GINI = mean(gini_disp_y0)), by = c("cntry", "essround")]
dim(table_A3)

write.csv(table_A3, file = "table_A3.csv")


##### Table A4: Welfare regime classification
table_A4_init <- table(as.character(ess_4_ld_models$cntry), ess_4_ld_models$welfare_regime)

table_A4 <- ifelse(test = table_A4_init == 0, yes = 0, no = 1)

write.csv(table_A4, file = "table_A4.csv")
